## Code to run fractional regression model for urban paper

## Started November 8, 2019.
## This version May 7th, 2021. Add in corrected county ids.

## clean up memory
rm(list=ls())

## Load in frm package
library(frm)

## Actually want frmhet
library(frmhet)

library(readstata13)

## load in data
## data <- read.csv("DataAugust132019.csv", h=T)
data <- read.dta13("Ready2020.dta")

data <- subset(data, (!is.na(us) & !is.na(spdnat) & !is.na(lnpopin1) &
						!is.na(dist) & !is.na(sharein1) & !is.na(protin1) & 
						!is.na(empl) & !is.na(alllandin1) &
                        !is.na(cathin1) & year==1949 & 
                        cathin1<100 & turnoutnat<100))

## Now basic FRM without accounting for endogeneity or fixed effects
## frm does not allow formula interface so have to specify y and x manually
y <- data$spdnat/100
X.names <- c("lnpopin1", "sharein1", "rel", "unempl", 
					   "alllandin1", "housing", "turnoutnat")
Z.names <- c("us", "dist", "dist2", "distus", "dist2us")
X <- data[, X.names]
Z <- data[, Z.names]

#####################################################
## Now estimate baseline model -- see if anything goes haywire.
## Logit
## Column (1) in Table B7
fr.model <- frm(y=y, x=X, linkfrac="logit")

AA <- frm.pe(fr.model)

##################################################
## Drop fixed effects
## This is Column (2) in Table B.7
Z.iv <- cbind(X[, -1], Z)

## Logit -- GMMz
fr.model.iv <- frmhet(y=y, x=X, z=Z.iv, type="GMMz", link="logit")

BB <- frmhet.pe(fr.model.iv)

## Now basic FRM without accounting for endogeneity or fixed effects
## frm does not allow formula interface so have to specify y and x manually
y <- data$cdunat/100
X.names <- c("lnpopin1", "sharein1", "rel", "unempl", 
					  "alllandin1", "housing", "turnoutnat")
Z.names <- c("us", "dist", "dist2", "distus", "dist2us")
X <- data[, X.names]
Z <- data[, Z.names]

####################################################
## Now estimate baseline model -- see if anything goes haywire.
## Logit
## This is Column (3) in Table B.7
fr.model <- frm(y=y, x=X, linkfrac="logit")

CC <- frm.pe(fr.model)

####################################################
## Drop fixed effects
## This is Column (4) in Table B7
X <- X[,1:2]
Z.iv <- cbind(Z, X[, -1])

## Logit -- GMMz
fr.model.iv <- frmhet(y=y, x=X, z=Z.iv, type="GMMz", link="cloglog")

DD <- frmhet.pe(fr.model.iv)

####################################################
## Construct Table B7

library(Hmisc)

table <- matrix(0,14,4)

sid1 <- seq(1, 14, by=2)
sid2 <- seq(2, 14, by=2)

table[sid1, 1] <- AA$PE.p[1:7]
table[sid2, 1] <- AA$PE.sd[1:7]

table[sid1, 2] <- BB$PE.p[1:7]
table[sid2, 2] <- BB$PE.sd[1:7]

table[sid1, 3] <- CC$PE.p[1:7]
table[sid2, 3] <- CC$PE.sd[1:7]

table[sid1, 4] <- DD$PE.p[1:7]
table[sid2, 4] <- DD$PE.sd[1:7]

name <- rep(" ", 14)
name[sid1] <- names(AA$PE.p)

rownames(table) <- name
table <- round(table, 4)

latex(table, file="table.tex")




















